/* "Efficiency and water use: Dynamic effects of irrigation technology adoption"
by Micah Cameron-Harp and Nathan Hendricks

Code written by Micah Cameron-Harp
May 16th, 2024

This do file creates the TWFE ATT and dynamic treatment effect estimates,
	and it also displays the sub-sample information in Table 2. Subsequent 
	do files combine these estimates with those produced using the 
	Callaway & Sant'Anna estimator in R to generate figures 3, 4, and 5.
*/

*Define directories and set working directory
	/* NOTE - To replicate our results, you need to change the directory address in the next line */
global dr_root = "C:\Users\Micah\Dropbox\Irrigation technology transition\final revisions for conditional acceptance\replication materials"
global dr_data = "$dr_root\data"
global dr_output = "$dr_root\outputs"
global dr_output_main = "$dr_root\outputs\main_text"
global dr_output_app = "$dr_root\outputs\appendices"
global dr_output_log = "$dr_root\outputs\logs"
global dr_temp = "$dr_root\data\intermediate"
cd "$dr_temp"


/* Start with Flood to CP or LEPA */
import excel "${dr_output_main}\twfe.xlsx", firstrow sheet(floodtocplepa_91_15_nyt) clear
	gen est = "twfe"
	replace scaled_estimate = 100*estimate/panel_mean_dep_var
	replace lb_estimate = estimate - 1.96*se_estimate
	replace ub_estimate = estimate + 1.96*se_estimate
	replace scaled_lb = 100*lb_estimate/panel_mean_dep_var
	replace scaled_ub = 100*ub_estimate/panel_mean_dep_var
	*Store the panel mean dep_var values for later
		sum panel_mean_dep_var if dep_var=="af_used"
		local mean_af_used = r(mean)
		sum panel_mean_dep_var if dep_var=="acres_irr"
		local mean_acres_irr = r(mean)
		sum panel_mean_dep_var if dep_var=="depth_applied"
		local mean_depth_applied = r(mean)
	keep dep_var est scaled_estimate scaled_lb scaled_ub estimate se_estimate lb_estimate ub_estimate
	save "${dr_temp}\floodtocplepa_twfe.dta", replace
	
import excel "${dr_temp}/cs_floodtocporlepa_91_15_nyt.xlsx", firstrow sheet(cs_att) clear
	keep if effect_estimated=="group_agg_effect_average"
	gen est = "cs"
	gen scaled_estimate = 100*estimate/panel_mean_dep_var
	gen scaled_lb = 100*ci_95_lb/panel_mean_dep_var
	gen scaled_ub = 100*ci_95_ub/panel_mean_dep_var
	rename (ci_95_lb ci_95_ub) (lb_estimate ub_estimate)
	keep dep_var est scaled_estimate scaled_lb scaled_ub estimate ///
		se_estimate lb_estimate ub_estimate
	save "${dr_temp}\floodtocplepa_cs.dta", replace
		
	*Generate column indicating tech transition and save
	gen transition = "flood_to_cp"
	save "${dr_output_tables}\floodtocplepa_att_allest.dta", replace 
	
use "${dr_temp}\didl_floodtocporlepa_afused.dta", clear
	gen dep_var = "af_used"
	gen panel_mean_dep_var=`mean_af_used'
	*Append the other two dependent variables 
	append using "${dr_output_tables}\flood_to_cplepa_acres_irr_didl_100breps.dta", force
	replace dep_var = "acres_irr" if missing(dep_var)
	replace panel_mean_dep_var=`mean_acres_irr' if missing(panel_mean_dep_var)
	append using "${dr_output_tables}\flood_to_cplepa_depth_applied_didl_100breps.dta", force
	replace dep_var = "depth_applied" if missing(dep_var)
	replace panel_mean_dep_var=`mean_depth_applied' if missing(panel_mean_dep_var)
	keep if time_to_treatment==.
	gen est = "cd"
	rename (treatment_effect se_treatment_effect treatment_effect_upper_95CI treatment_effect_lower_95CI) ///
		(estimate se_estimate ub_estimate lb_estimate)
	gen scaled_estimate = 100*estimate/panel_mean_dep_var
	gen scaled_lb = 100*lb_estimate/panel_mean_dep_var
	gen scaled_ub = 100*ub_estimate/panel_mean_dep_var
	keep dep_var est scaled_estimate scaled_lb scaled_ub estimate se_estimate lb_estimate ub_estimate

	*Append the TWFE and CS results 
	append using "${dr_temp}\floodtocplepa_twfe.dta"
	append using "${dr_temp}\floodtocplepa_cs.dta"
	*Drop "af_used_irr" dep_var
	drop if dep_var=="af_used_irr"	
		
	*Generate column indicating tech transition and save
	gen transition = "flood_to_cp"
	save "${dr_output_tables}\floodtocplepa_att_allest.dta", replace 

	
/* Now CP to LEPA */
import excel "${dr_output_tables}\twfe.xlsx", firstrow sheet(cp to lepa) clear
	gen est = "twfe"
	replace scaled_estimate = 100*estimate/panel_mean_dep_var
	replace lb_estimate = estimate - 1.96*se_estimate
	replace ub_estimate = estimate + 1.96*se_estimate
	gen var_estimate = se_estimate^2
	gen scaled_se = (var_estimate*(100/panel_mean_dep_var)^2)^.5
	replace scaled_lb = scaled_estimate - 1.96*scaled_se
	replace scaled_ub = scaled_estimate + 1.96*scaled_se
	*Store the panel mean dep_var values for later
		sum panel_mean_dep_var if dep_var=="af_used"
		local mean_af_used = r(mean)
		sum panel_mean_dep_var if dep_var=="acres_irr"
		local mean_acres_irr = r(mean)
		sum panel_mean_dep_var if dep_var=="depth_applied"
		local mean_depth_applied = r(mean)
	keep dep_var est scaled_estimate scaled_lb scaled_ub estimate se_estimate lb_estimate ub_estimate
	save "${dr_temp}\cptolepa_twfe.dta", replace
	
import excel "${r_dr_output_tables}/cs_cptolepa.xlsx", firstrow sheet(cs_att) clear
	keep if effect_estimated=="group_agg_effect_average"
	gen est = "cs"
	gen scaled_estimate = 100*estimate/panel_mean_dep_var
	gen scaled_lb = 100*ci_95_lb/panel_mean_dep_var
	gen scaled_ub = 100*ci_95_ub/panel_mean_dep_var
	rename (ci_95_lb ci_95_ub) (lb_estimate ub_estimate)
	keep dep_var est scaled_estimate scaled_lb scaled_ub estimate se_estimate lb_estimate ub_estimate
	save "${dr_temp}\cptolepa_cs.dta", replace
	
use "${dr_output_tables}\cp_to_lepa_af_used_didl_100breps.dta", clear
	gen dep_var = "af_used"
	gen panel_mean_dep_var=`mean_af_used'
	*Append the other two dependent variables 
	append using "${dr_output_tables}\cp_to_lepa_acres_irr_didl_100breps.dta", force
	replace dep_var = "acres_irr" if missing(dep_var)
	replace panel_mean_dep_var=`mean_acres_irr' if missing(panel_mean_dep_var)
	append using "${dr_output_tables}\cp_to_lepa_depth_applied_didl_100breps.dta", force
	replace dep_var = "depth_applied" if missing(dep_var)
	replace panel_mean_dep_var=`mean_depth_applied' if missing(panel_mean_dep_var)
	keep if time_to_treatment==.
	gen est = "cd"
	rename (treatment_effect se_treatment_effect treatment_effect_upper_95CI treatment_effect_lower_95CI) ///
		(estimate se_estimate ub_estimate lb_estimate)
	gen scaled_estimate = 100*estimate/panel_mean_dep_var
	gen scaled_lb = 100*lb_estimate/panel_mean_dep_var
	gen scaled_ub = 100*ub_estimate/panel_mean_dep_var
	keep dep_var est scaled_estimate scaled_lb scaled_ub estimate se_estimate lb_estimate ub_estimate

	*Append the TWFE and CS results 
	append using "${dr_temp}\cptolepa_twfe.dta"
	append using "${dr_temp}\cptolepa_cs.dta"
	*Drop "af_used_irr" dep_var
	drop if dep_var=="af_used_irr"	
		
	*Generate column indicating tech transition and save
	gen transition = "cp_to_lepa"
	save "${dr_output_tables}\cptolepa_att_allest.dta", replace 

********************************************************************************
*Combine all estimated ATTs for all estimators and transitions
	append using "${dr_output_tables}\floodtocplepa_att_cs_twfe.dta"
	*Make numerical variable encoding estimator 
	gen estimator = 1 if est=="cs"
	replace estimator = 2 if est=="twfe"
	drop if est=="cd"

/************************************* PLOTTING *******************************/

	*Acre-ft, estimates that passed placebo in navy, failed in red
	tw (rcap scaled_lb scaled_ub estimator if transition=="flood_to_cp" & ///
		dep_var=="af_used", lcolor("navy")) ///
		(scatter scaled_estimate estimator if transition=="flood_to_cp" & ///
		dep_var=="af_used", mcolor("navy") msize(small)), ///  
		title("Acre-feet withdrawn", size(mediumsmall)) ///
		yscale(range(-15 15)) ylabel(-15(5)15, labsize(small)) xscale(range(.8 2.2)) ///
		legend(off) ///
		yline(0, lpattern(solid)) ///
		ytitle("") xtitle("") ///
		xlabel(1 "{it}{&delta}{superscript:CS}" 2 "{it}{&beta}{superscript:TWFE}", labsize(medium)) ///
		graphregion(color(white) margin (1 1 1 8)) ///
		text(23 .3 "A. Flood to center pivot", place(e) size(medium)) ///
		xsize(2.5) ysize(3) ///
		name(flood_cp_acreft, replace)		
		
		// 		legend(order(2) lab(2 "Estimate as percent of the dependent variable's mean value")) ///
		
	*Acres irrigated, estimates that passed placebo in navy, failed in red
	tw (rcap scaled_lb scaled_ub estimator if transition=="flood_to_cp" & ///
		dep_var=="acres_irr", lcolor("navy")) ///
		(scatter scaled_estimate estimator if transition=="flood_to_cp" & ///
		dep_var=="acres_irr", mcolor("navy") msize(small)), ///  
		title("Acres irrigated", size(mediumsmall)) ///
		yscale(range(-15 15)) ylabel(-15(5)15, labsize(small)) xscale(range(.8 2.2)) ///
		legend(off) ///
		yline(0, lpattern(solid)) ///
		ytitle("") xtitle("") ///
		xlabel(1 "{it}{&delta}{superscript:CS}" 2 "{it}{&beta}{superscript:TWFE}", labsize(medium)) ///
		graphregion(color(white) margin (1 1 1 8)) ///
		xsize(2) ysize(3) ///
		name(flood_cp_acresirr, replace)
		
	*Depth applied, estimates that passed placebo in navy, failed in red
	tw (rcap scaled_lb scaled_ub estimator if transition=="flood_to_cp" & ///
		dep_var=="depth_applied", lcolor("navy")) ///
		(scatter scaled_estimate estimator if transition=="flood_to_cp" & ///
		dep_var=="depth_applied", mcolor("navy") msize(small)), ///  
		title("Depth applied", size(mediumsmall)) ///
		yscale(range(-15 15)) ylabel(-15(5)15, labsize(small)) xscale(range(.8 2.2)) ///
		legend(off) ///
		yline(0, lpattern(solid)) ///
		ytitle("") xtitle("") ///
		xlabel(1 "{it}{&delta}{superscript:CS}" 2 "{it}{&beta}{superscript:TWFE}", labsize(medium)) ///
		graphregion(color(white) margin (1 1 1 8)) ///
		xsize(2) ysize(3) ///
		name(flood_cp_depthapp, replace)
		
	*Combine all the flood to cp graphs 
		graph combine flood_cp_acreft flood_cp_acresirr flood_cp_depthapp, ///
		rows(1) cols(3) ///
		l1title("Percent change", size(small)) ///
		xsize(6.5) ysize(2.65) name(flood_cp_all, replace)
			
**# CP to cp with Lepa
	*Acre-ft, estimates that passed placebo in navy, failed in red
	tw (rcap scaled_lb scaled_ub estimator if transition=="cp_to_lepa" & ///
		dep_var=="af_used", lcolor("navy")) ///
		(scatter scaled_estimate estimator if transition=="cp_to_lepa" & ///
		dep_var=="af_used", mcolor("navy") msize(small)), ///  
		title("Acre-feet withdrawn", size(mediumsmall)) ///
		yscale(range(-10 5)) ylabel(-10(2.5)5, labsize(small)) xscale(range(.8 2.2)) ///
		legend(off) ///
		yline(0, lpattern(solid)) ///
		ytitle("") xtitle("") ///
		xlabel(1 "{it}{&delta}{superscript:CS}" 2 "{it}{&beta}{superscript:TWFE}", labsize(medium)) ///
		graphregion(color(white) margin (1 1 1 8)) ///
		text(8.8 .29 "B. Traditional center pivot to LEPA", place(e) size(medium)) ///
		xsize(2.5) ysize(3) ///
		name(cp_lepa_acreft, replace) 
		
	*Acres irrigated, estimates that passed placebo in navy, failed in red
	tw (rcap scaled_lb scaled_ub estimator if transition=="cp_to_lepa" & ///
		dep_var=="acres_irr", lcolor("navy")) ///
		(scatter scaled_estimate estimator if transition=="cp_to_lepa" & ///
		dep_var=="acres_irr", mcolor("navy") msize(small)), ///  
		title("Acres irrigated", size(mediumsmall)) ///
		yscale(range(-10 5)) ylabel(-10(2.5)5, labsize(small)) xscale(range(.8 2.2)) ///
		legend(off) ///
		yline(0, lpattern(solid)) ///
		ytitle("") xtitle("") ///
		xlabel(1 "{it}{&delta}{superscript:CS}" 2 "{it}{&beta}{superscript:TWFE}", labsize(medium)) ///
		graphregion(color(white) margin (1 1 1 8)) /// 
		xsize(2) ysize(3) ///
		name(cp_lepa_acresirr, replace)
		
	*Depth applied, estimates that passed placebo in navy, failed in red
	tw (rcap scaled_lb scaled_ub estimator if transition=="cp_to_lepa" & ///
		dep_var=="depth_applied", lcolor("navy")) ///
		(scatter scaled_estimate estimator if transition=="cp_to_lepa" & ///
		dep_var=="depth_applied", mcolor("navy") msize(small)), ///  
		title("Depth applied", size(mediumsmall)) ///
		yscale(range(-10 5)) ylabel(-10(2.5)5, labsize(small)) xscale(range(.8 2.2)) ///
		legend(off) ///
		yline(0, lpattern(solid)) ///
		ytitle("") xtitle("") ///
		xlabel(1 "{it}{&delta}{superscript:CS}" 2 "{it}{&beta}{superscript:TWFE}", labsize(medium)) ///
		graphregion(color(white) margin (1 1 1 8)) ///
		xsize(2) ysize(3) ///
		name(cp_lepa_depthapp, replace)
		
	*Combine all the cp to LEPA graphs
		graph combine cp_lepa_acreft cp_lepa_acresirr cp_lepa_depthapp, ///
		rows(1) cols(3) ///
		l1title("Percent change", size(small)) ///
		xsize(6.5) ysize(2.65) name(cp_lepa_all, replace)
	
*Combine the graphs 	
graph combine flood_cp_all cp_lepa_all, ///
	rows(2) cols(1) ///
	b1title("Estimator", size(small)) ///
	xsize(6.5) ysize(8) 
*Open graph editor and move row labels
graph export "C:\Users\Micah\Dropbox\Irrigation technology transition\2nd round of revisions\final code\stata_do\output\graphs\ONLY CS and TWFE_att_scaled estimates_final version.tif", width(6500) height(7800) replace 

